Hygrogen Atom Problem - Analytical Solutions and Graphs - SKP¶

In [1]:
import numpy as np
import scipy as scp
import matplotlib.pyplot as plt
plt.style.use(['dark_background'])
import sympy as sp
import sympy.physics.hydrogen as sphydrogen
In [2]:
import scipy.constants as const

$E_h = 27.21138602 \,eV\, = 4.35974434 × 10^{-18} \, J$

Description¶

  • I have defined the functions in the Defining Function section. Here I have taken functions from sympy.physics.hydrogen and created numerical functions from them. The defined function, named hydrogen_wavefn_all gives the analytical and graphical view of total wavefunction $\psi_{nlm}(r, \theta, \phi)$ and radial wavefunction $R_{nl}(r)$ and hydrogen_wavefn_return returns 2 functions phi_nlm_fn(r, theta, phi) and R_nl_fn(r) for those wavefunctions.

  • In the Give Inputs and Get States section, you can give arbitary values of (n, l, m) and get the states. Here you need to set the number of grid points (N1) and the maximum limit of radius r1 to view the total graphs.

  • Once you are done with the above things, go to the Working with Functions section. Here I have used the returned functions (phi_nlm_fn(r, theta, phi) and R_nl_fn(r)) and done few more plots from these.

Defining Function¶

Function: hydrogen_wavefn_all¶

cmap: 'Accent', 'Accent_r', 'Blues', 'Blues_r', 'BrBG', 'BrBG_r', 'BuGn', 'BuGn_r', 'BuPu', 'BuPu_r', 'CMRmap', 'CMRmap_r', 'Dark2', 'Dark2_r', 'GnBu', 'GnBu_r', 'Greens', 'Greens_r', 'Greys', 'Greys_r', 'OrRd', 'OrRd_r', 'Oranges', 'Oranges_r', 'PRGn', 'PRGn_r', 'Paired', 'Paired_r', 'Pastel1', 'Pastel1_r', 'Pastel2', 'Pastel2_r', 'PiYG', 'PiYG_r', 'PuBu', 'PuBuGn', 'PuBuGn_r', 'PuBu_r', 'PuOr', 'PuOr_r', 'PuRd', 'PuRd_r', 'Purples', 'Purples_r', 'RdBu', 'RdBu_r', 'RdGy', 'RdGy_r', 'RdPu', 'RdPu_r', 'RdYlBu', 'RdYlBu_r', 'RdYlGn', 'RdYlGn_r', 'Reds', 'Reds_r', 'Set1', 'Set1_r', 'Set2', 'Set2_r', 'Set3', 'Set3_r', 'Spectral', 'Spectral_r', 'Wistia', 'Wistia_r', 'YlGn', 'YlGnBu', 'YlGnBu_r', 'YlGn_r', 'YlOrBr', 'YlOrBr_r', 'YlOrRd', 'YlOrRd_r', 'afmhot', 'afmhot_r', 'autumn', 'autumn_r', 'binary', 'binary_r', 'bone', 'bone_r', 'brg', 'brg_r', 'bwr', 'bwr_r', 'cividis', 'cividis_r', 'cool', 'cool_r', 'coolwarm', 'coolwarm_r', 'copper', 'copper_r', 'cubehelix', 'cubehelix_r', 'flag', 'flag_r', 'gist_earth', 'gist_earth_r', 'gist_gray', 'gist_gray_r', 'gist_heat', 'gist_heat_r', 'gist_ncar', 'gist_ncar_r', 'gist_rainbow', 'gist_rainbow_r', 'gist_stern', 'gist_stern_r', 'gist_yarg', 'gist_yarg_r', 'gnuplot', 'gnuplot2', 'gnuplot2_r', 'gnuplot_r', 'gray', 'gray_r', 'hot', 'hot_r', 'hsv', 'hsv_r', 'inferno', 'inferno_r', 'jet', 'jet_r', 'magma', 'magma_r', 'nipy_spectral', 'nipy_spectral_r', 'ocean', 'ocean_r', 'pink', 'pink_r', 'plasma', 'plasma_r', 'prism', 'prism_r', 'rainbow', 'rainbow_r', 'seismic', 'seismic_r', 'spring', 'spring_r', 'summer', 'summer_r', 'tab10', 'tab10_r', 'tab20', 'tab20_r', 'tab20b', 'tab20b_r', 'tab20c', 'tab20c_r', 'terrain', 'terrain_r', 'turbo', 'turbo_r', 'twilight', 'twilight_r', 'twilight_shifted', 'twilight_shifted_r', 'viridis', 'viridis_r', 'winter', 'winter_r'

Selected: 'viridis', 'plasma', 'inferno', 'magma'

In [3]:
def hydrogen_wavefn_all(Z, n, l, m, r, theta, phi):
    '''
    ------- Functions to import --------
    import numpy as np
    import matplotlib.pyplot as plt
    import sympy as sp
    import sympy.physics.hydrogen as sphydrogen
    -------- Inputs ---------
    Z: atomic number
    n: principal quantum number: 1,2,3,...
    l: azimuthal quantum number: 0,1,2,... (l<n-1)
    m: magnetic quantum number: ...,-2,-1,0,1,2,... (|m|<l)
    r: an numpy array: np.linspace(0,15,N)
    theta: numpy array: np.linspace(0, np.pi, N)
    phi: numpy array: np.linspace(0, 2*np.pi, N1)
    N = 50  # number of grid points for plotting
    ------- Outputs ----------
    value of energy of the state
    Analytical expression for psi_nlm(r,theta,phi)
    Normalization check of psi_nlm(r, theta, phi) (compare with 1)
    Contour plot for psi in (r*cos(psi), r*sin(phi)) plane
    Analytical expression for R_nl(r)
    Normalization check of R_nl(r) (compare with 1)
    Plot of r vs R_nl(r)
    '''
    print(f'Inputs: n={n}, l={l}, m={m}, Z={Z}')
    rsp = sp.Symbol('r', positive=True)
    phisp, thetasp = sp.symbols('phi theta', real=True)
    Zsp = sp.Symbol('Z', positive=True, integer=True, nonzero=True)
    psinlmsp, rnlsp = sp.symbols('\psi_{nlm} R_{nl}')

    EheV, EhJ = 27.21138602, 4.35974434e-18  # Energy in eV and J
    E_n = sphydrogen.E_nl(n, Z)
    print(f'E = {E_n} in Hartree atomic units = {EhJ*E_n} J = {EheV*E_n} eV')

    Psi_nlm_sp = sphydrogen.Psi_nlm(n, l, m, rsp, phisp, thetasp, Zsp)
    display(psinlmsp, Psi_nlm_sp, Psi_nlm_sp.simplify())
    print('normalization:', sp.integrate(sp.Abs(Psi_nlm_sp)**2*rsp**2*sp.sin(thetasp), 
                    (rsp,0,sp.oo),(thetasp,0,sp.pi),(phisp,0,2*sp.pi)))
    Psi_nlm_np = sp.lambdify((rsp,phisp,thetasp), sphydrogen.Psi_nlm(n,l,m,rsp,phisp,thetasp,Z))
    # R3, Theta3, Phi3 = np.meshgrid(r, theta, phi, indexing='ij')
    # psiarr3 = Psi_nlm_np(R3, Theta3, Phi3)
    R2, Phi2 = np.meshgrid(r, phi)
    Theta2 = np.full_like(R2, np.pi/2)  # input value of theta
    psiarr2 = Psi_nlm_np(R2, Theta2, Phi2)
    plt.contourf(R2*np.cos(Phi2), R2*np.sin(Phi2), psiarr2, cmap='inferno')
    plt.colorbar(label='$\psi_{nlm}$')
    plt.title(f'State {n, l, m}, Z={Z}')
    plt.show()
    plt.contourf(R2*np.cos(Phi2), R2*np.sin(Phi2), np.abs(psiarr2)**2, cmap='magma')
    plt.colorbar(label='$|\psi_{nlm}|^2$')
    plt.title(f'State {n, l, m}, Z={Z}')
    plt.show()

    R_nl_sp = sphydrogen.R_nl(n, l, rsp, Zsp)
    display(rnlsp, R_nl_sp, R_nl_sp.simplify())
    R_nl_np = sp.lambdify(rsp, sphydrogen.R_nl(n, l, rsp, Z))
    R_nl_plot = R_nl_np(r)/(np.sum(R_nl_np(r)**2*r**2*(r[1]-r[0])))**0.5
    print(f'normalization: {np.sum(R_nl_plot**2 *r**2 *(r[1]-r[0]))}')
    plt.plot(r, R_nl_plot)
    plt.axhline(0, color='yellow')
    plt.xlabel('$r$')
    plt.ylabel('$R_{nl}$', fontsize=12)
    plt.title(f'State {n, l}, Z={Z}')
    plt.grid()
    plt.show()
    plt.plot(r, R_nl_plot**2)
    plt.axhline(0, color='yellow')
    plt.xlabel('$r$')
    plt.ylabel('$|R_{nl}|^2$', fontsize=12)
    plt.title(f'State {n, l}, Z={Z}')
    plt.grid()
    plt.show()
In [4]:
print(hydrogen_wavefn_all.__doc__)

    ------- Functions to import --------
    import numpy as np
    import matplotlib.pyplot as plt
    import sympy as sp
    import sympy.physics.hydrogen as sphydrogen
    -------- Inputs ---------
    Z: atomic number
    n: principal quantum number: 1,2,3,...
    l: azimuthal quantum number: 0,1,2,... (l<n-1)
    m: magnetic quantum number: ...,-2,-1,0,1,2,... (|m|<l)
    r: an numpy array: np.linspace(0,15,N)
    theta: numpy array: np.linspace(0, np.pi, N)
    phi: numpy array: np.linspace(0, 2*np.pi, N1)
    N = 50  # number of grid points for plotting
    ------- Outputs ----------
    value of energy of the state
    Analytical expression for psi_nlm(r,theta,phi)
    Normalization check of psi_nlm(r, theta, phi) (compare with 1)
    Contour plot for psi in (r*cos(psi), r*sin(phi)) plane
    Analytical expression for R_nl(r)
    Normalization check of R_nl(r) (compare with 1)
    Plot of r vs R_nl(r)

Function: hydrogen_wavefn_return¶

In [5]:
def hydrogen_wavefn_return(Z, n, l, m):
    '''
    ------- Functions to import --------
    import numpy as np
    import sympy as sp
    import sympy.physics.hydrogen as sphydrogen
    -------- Inputs ---------
    Z: atomic number
    n: principal quantum number: 1,2,3,...
    l: azimuthal quantum number: 0,1,2,... (l<n-1)
    m: magnetic quantum number: ...,-2,-1,0,1,2,... (|m|<l)
    ------- Returns ----------
    np.array([E, E in J, E in eV]) psi_nlm(r,theta,phi), R_nl(r)
    '''
    # print(f'Inputs: n={n}, l={l}, m={m}, Z={Z}')
    rsp = sp.Symbol('r', positive=True)
    phisp, thetasp = sp.symbols('phi theta', real=True)

    EheV, EhJ = 27.21138602, 4.35974434e-18  # Energy in eV and J
    E_n = sphydrogen.E_nl(n, Z)
    E_n_J, E_n_eV = EhJ*E_n, EheV*E_n
    Psi_nlm_np = sp.lambdify((rsp,phisp,thetasp), sphydrogen.Psi_nlm(n,l,m,rsp,phisp,thetasp,Z))
    R_nl_np = sp.lambdify(rsp, sphydrogen.R_nl(n, l, rsp, Z))

    return np.array([E_n, E_n_J, E_n_eV]), Psi_nlm_np, R_nl_np
In [6]:
print(hydrogen_wavefn_return.__doc__)

    ------- Functions to import --------
    import numpy as np
    import sympy as sp
    import sympy.physics.hydrogen as sphydrogen
    -------- Inputs ---------
    Z: atomic number
    n: principal quantum number: 1,2,3,...
    l: azimuthal quantum number: 0,1,2,... (l<n-1)
    m: magnetic quantum number: ...,-2,-1,0,1,2,... (|m|<l)
    ------- Returns ----------
    np.array([E, E in J, E in eV]) psi_nlm(r,theta,phi), R_nl(r)

Give Inputs and Get States¶

In [7]:
N1 = 500
rmax = 50
Z1 = 1
n1, l1, m1 = 5,3,-1

r1 = np.linspace(0, rmax, N1)
theta1 = np.linspace(0, np.pi, N1)
phi1 = np.linspace(0, 2*np.pi, N1)
hydrogen_wavefn_all(Z1, n1,l1,m1,r1,theta1,phi1)

Inputs: n=5, l=3, m=-1, Z=1
E = -1/50 in Hartree atomic units = -8.71948868000000E-20 J = -0.544227720400000 eV
$\displaystyle \psi_{nlm}$
$\displaystyle \frac{\sqrt{70} Z^{\frac{9}{2}} r^{3} \left(- \frac{2 Z r}{5} + 8\right) \left(\frac{5 \sqrt{21} e^{- i \phi} \sin{\left(\theta \right)} \cos^{2}{\left(\theta \right)}}{8 \sqrt{\pi}} - \frac{\sqrt{21} e^{- i \phi} \sin{\left(\theta \right)}}{8 \sqrt{\pi}}\right) e^{- \frac{Z r}{5}}}{328125}$
$\displaystyle - \frac{\sqrt{30} Z^{\frac{9}{2}} r^{3} \left(Z r - 20\right) \left(5 \cos^{2}{\left(\theta \right)} - 1\right) e^{- \frac{Z r}{5} - i \phi} \sin{\left(\theta \right)}}{937500 \sqrt{\pi}}$ 
normalization: 1

C:\ProgramData\Anaconda3\lib\site-packages\numpy\ma\core.py:2829: ComplexWarning: Casting complex values to real discards the imaginary part
  _data = np.array(data, dtype=dtype, copy=copy,
$\displaystyle R_{nl}$
$\displaystyle \frac{\sqrt{70} Z^{\frac{9}{2}} r^{3} \left(- \frac{2 Z r}{5} + 8\right) e^{- \frac{Z r}{5}}}{328125}$
$\displaystyle \frac{2 \sqrt{70} Z^{\frac{9}{2}} r^{3} \left(- Z r + 20\right) e^{- \frac{Z r}{5}}}{1640625}$ 
normalization: 1.0

Working with Functions¶

Main Function: hydrogen_wavefn_return(Z, n,l,m)

Returns:

  1. $E_n \, \implies \,$ np.array([E, E in J, E in eV])
  2. $\psi_{nlm}(r, \theta, \phi) \, \implies \,$ phi_nlm_fn(r, theta, phi)
  3. $R_{nl}(r) \, \implies \,$ R_nl_fn(r)
In [8]:
Z1, n1, l1, m1 = 1, 2, 1, -1   # INPUT
En, psi_nlm_fn, R_nl_fn = hydrogen_wavefn_return(Z1, n1,l1,m1)
display(En, psi_nlm_fn, R_nl_fn)

array([-1/8, -5.44968042500000e-19, -3.40142325250000], dtype=object)
<function _lambdifygenerated(r, phi, theta)>
<function _lambdifygenerated(r)>

$\psi_{nlm}(r, \theta, \phi)$ in plotly¶

In [9]:
import plotly.graph_objects as go
In [10]:
N3 = 200
r3 = np.linspace(0, n1*12, N3)  # input
phi3 = np.linspace(0, 2*np.pi, N3)
R2, Phi2 = np.meshgrid(r3, phi3)
theta_const = 0.5  # input
Theta2 = np.full_like(R2, theta_const)
psinlmplotly = psi_nlm_fn(R2, Theta2, Phi2)

fig1 = go.Figure(data=[go.Surface(z=np.abs(psinlmplotly)**2, 
            x=R2*np.cos(Phi2), y=R2*np.sin(Phi2), colorscale='turbo')])
fig1.update_layout(scene=dict(
                xaxis_title='r cos(phi)',
                yaxis_title='r sin(phi)',
                zaxis_title='|psi_nlm|^2'),
                title=f'State {n1,l1,m1}')
fig1.show()
In [11]:
N4 = 20
r4 = np.linspace(0, n1*12, N4)  # input
theta4 = np.linspace(0, np.pi, N4)
phi4 = np.linspace(0, 2*np.pi, N4)
R3, Theta3, Phi3 = np.meshgrid(r4, theta4, phi4)
psi2nlmplotly = np.abs(psi_nlm_fn(R3, Theta3, Phi3))**2
X3, Y3 = R3*np.sin(Theta3)*np.cos(Phi3), R3*np.sin(Theta3)*np.sin(Phi3)
Z3 = R3*np.cos(Theta3)

isovalue = np.abs(psi2nlmplotly).max() *0.5
fig2 = go.Figure(data=go.Isosurface(
    x = X3.flatten(),
    y = Y3.flatten(),
    z = Z3.flatten(),
    value = psi2nlmplotly.flatten(),
    isomin=0,
    isomax=isovalue,
    opacity = 0.7, # CHANGE
    surface_count=2,
    colorscale='plasma' ))

fig2.update_layout(
                title=f'State {n1,l1,m1}, Energy={En[2]:.4f}eV',
                scene=dict(xaxis_title='x', yaxis_title='y', zaxis_title='z'))
fig2.show()
In [12]:
display(psi2nlmplotly.flatten().shape, Z3.flatten().shape)

(8000,)
(8000,)

All $R_{nl}(r)$ plots for a particular $n$¶

In [13]:
Z2, n2 = 1, 5
rmax2, N2 = (n2+2)*10, 1000
r2 = np.linspace(0, rmax2, N2)
Rnl_all = list(np.zeros(n2))
for l2 in range(n2):
    En2, psinlm2, Rnl_all[l2] = hydrogen_wavefn_return(Z2,n2,l2,0)
    Rnl_plot = np.array(Rnl_all[l2](r2))
    print(f'normalization: {np.sum(Rnl_plot**2*r2**2*(r2[1]-r2[0]))}')
    plt.plot(r2, Rnl_plot)
    plt.axhline(0, color='yellow')
    plt.xlabel('$r$')
    plt.ylabel('$R_{nl}$', fontsize=12)
    plt.title(f'State {n2, l2}, Z={Z2}')
    plt.grid()
    plt.show()

normalization: 0.9978348220207655

normalization: 0.9983043342534061

normalization: 0.9989981578914803

normalization: 0.9995956114985672

normalization: 0.999915281662709
In [ ]: